#ifndef LB_Sparse_Matrix_H_
#define LB_Sparse_Matrix_H_
#define Sparse_Matrix LB_Sparse_Matrix 
#include<Math/LB_Math.h>

#ifdef __cplusplus
extern "C" {
#endif

typedef struct SparseMatrixEntry
{
    int i;double v;
}SparseMatrixEntry;

static inline void sparsematrixentry_init(SparseMatrixEntry*sme)
{
    sme->i=0;sme->v=0;
}

typedef struct LB_Sparse_Matrix{
    int rows,cols;
    int* rowSizes;
    SparseMatrixEntry** data;

    void *prop;
}LB_Sparse_Matrix;

static inline void lb_sparse_matrix_initn(LB_Sparse_Matrix*lsm,int rows,int cols)
{
    lsm->rows=rows;lsm->cols=cols;
    lsm->rowSizes=(int*)malloc(sizeof(int)*rows);
    memset(lsm->rowSizes,0,sizeof(int)*rows);
    lsm->data=(SparseMatrixEntry**)malloc(sizeof(SparseMatrixEntry*)*rows);
    memset(lsm->data,0,sizeof(SparseMatrixEntry*)*rows);
    lsm->prop=NULL;
}

void lb_sparse_matrix_free(LB_Sparse_Matrix*lsm);


void lb_sparse_matrix_resize(LB_Sparse_Matrix*lsm,int rows,int cols);


static inline void   lb_sparse_matrix_setrowsize(LB_Sparse_Matrix*lsm,int row,int count)
{
    SAFE_FREE(lsm->data[row]);
    lsm->data[row]=(SparseMatrixEntry*)malloc(sizeof(SparseMatrixEntry)*count);
    lsm->rowSizes[row]=count;
}
LB_Sparse_Matrix* lb_sparse_matrix_transpose(LB_Sparse_Matrix* lsm);

// lsm 可以是re
void lb_sparse_matrix_multply_scala(LB_Sparse_Matrix* lsm,
    double fac,LB_Sparse_Matrix*re);

// copy lsm2 to lsm1
void lb_sparse_matrix_copy_data(LB_Sparse_Matrix* lsm1,LB_Sparse_Matrix* lsm2);



LB_Sparse_Matrix* lb_sparse_matrix_new_copy(LB_Sparse_Matrix* lsm);

// result 不能是lsm1或者lsm2
// 
void lb_sparse_matrix_multiply(LB_Sparse_Matrix*result,LB_Sparse_Matrix* lsm1, 
    LB_Sparse_Matrix* lsm2);
// result 可以等于lsm1 或者lsm2

void lb_sparse_matrix_plus(LB_Sparse_Matrix* result,
    LB_Sparse_Matrix* lsm1, 
    LB_Sparse_Matrix* lsm2);


// result 可以等于lsm1 或者lsm2
// 
void lb_sparse_matrix_minus(LB_Sparse_Matrix* result,
    LB_Sparse_Matrix* lsm1, 
    LB_Sparse_Matrix* lsm2);

void lb_sparse_matrix_mult_array(LB_Sparse_Matrix* lsm, double* vec,double* result); 


void lb_sparse_matrix_printself(LB_Sparse_Matrix *lsm);

// 
// 要保证lsm是对陈正定矩阵
//
void compute_conjugate_gradient(LB_Sparse_Matrix*lsm,double* b,double * x ,int iter,double eps);




// Sparse LU分解
void compute_lu_decompose(LB_Sparse_Matrix*lsm,
    LB_Sparse_Matrix* P,LB_Sparse_Matrix* L,LB_Sparse_Matrix*U);




// // Solve for x s.t. M(x)=b by solving for x s.t. M^tM(x)=M^t(b)
// template<class T>
// int SparseMatrix<T>::Solve(const SparseMatrix<T>& M,const Vector<T>& b,const int& iters,Vector<T>& solution,const T eps){
//     SparseMatrix mTranspose=M.Transpose();
//     Vector<T> bb=mTranspose*b;
//     Vector<T> d,r,Md;
//     T alpha,beta,rDotR;
//     int i;

//     solution.Resize(M.Columns());
//     solution.SetZero();

//     d=r=bb;
//     rDotR=r.Dot(r);

//     for(i=0;i<iters && rDotR>eps;i++){
//         T temp;
//         Md=mTranspose*(M*d);
//         alpha=rDotR/d.Dot(Md);
//         solution+=d*alpha;
//         r-=Md*alpha;
//         temp=r.Dot(r);
//         beta=temp/rDotR;
//         rDotR=temp;
//         d=r+d*beta;
//     }
//     return i;
// }



#ifdef __cplusplus
}
#endif

#undef Sparse_Matrix
#endif